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Abstract. Coherent-state quantum process tomography (csQPT) is a method of 
completely characterizing a quantum-optical "black box" by probing it with coherent 
states and performing homodyne measurements on the output [M. Lobino et al , 
Science 322, 563 (2008)]. We present a technique for csQPT that is fully based 
on statistical inference, specifically, quantum expectation-maximization. The method 
relies on the Jamiolkowski isomorphism and iteratively reconstructs the process tensor 
in the Fock basis directly from the experimental data. This approach permits 
incorporation of a priori constraints into the reconstruction procedure, thereby 
guaranteeing that the resulting process tensor is physically consistent. Furthermore, 
our method is easier to implement and requires a narrower range of coherent states 
than its predecessors. We test its feasibility using simulations on several experimentally 
relevant processes. 
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1. Introduction 

The art of determining states of quantum systems — quantum tomography — relies on 
performing measurements over multiple copies of the state in various bases, followed by 
reconstruction of the state's density matrix using suitable algorithms on the procured 
data. Methods of state tomography can be extended to the quantum version of the 
"black box" problem [TJ [2j [3], giving rise to quantum process tomography (QPT). In 
QPT, measurements on the black box response to a certain set of probe states allow 
one to predict the effect of that black box on any arbitrary state within a given Hilbert 
space. QPT emerged in response to ever-increasing demands in the field of quantum 
information processing, as the assembly of any quantum information processor requires 
precise knowledge of each of its components [I] . 

A popular approach to QPT involves determining the output £{pi) for each state 
of a spanning set {pi} of the space of density matrices over the Hilbert space of interest. 
Due to the linearity of quantum processes over its density operators, the output of any 
arbitrary state p = qpj can then be found as S(p) = J2i c i£(f>i)- 

This approach has recently been extended to the continuous-variable domain of 
quantum optics [5]. The reconstruction procedure involves probing the process with 
coherent states, i.e. simple laser pulses. It relies on the ability of the Glauber-Sudarshan 
P representation to express the density matrix of any quantum state as a linear 
combination of coherent states' density matrices. Improvements in the algorithm have 
been presented in [6] . The algorithm has been tested in an experiment on characterizing 
quantum-optical memory [7]. Similar principles have recently been used to perform 
characterization of quantum optical detectors [HJ [9] . 

This method, known as coherent-state QPT, or csQPT, has the advantage of 
employing only the easy-to-prepare coherent states for probing. However, the numerical 
reconstruction procedures employed in [5J E] involve an intermediate step of determining 
the density matrices of the output states £(|a)(a:|) for each probe coherent state |a) 
and subsequent integration with the P function. This approach requires a multistep 
calculation and does not guarantee to yield a process that is physically plausible, i.e. 
completely positive and trace non-increasing. 

We present a reconstruction scheme that does away with this intermediate step, 
and reconstructs the process directly from the experimental data using pure statistical 
inference. The experimental setup is equivalent to that of [5] and is illustrated in figure HJ 
The process reconstruction algorithm, on the other hand, is entirely different: it relies 
on the iterative maximum-likelihood approach. Its major advantage is the possibility 
to incorporate a priori constraints in the reconstruction procedure in order to ensure 
physically consistent and meaningful results. 

Maximum-likelihood methods have been successfully used in the past for quantum 
state estimation as well as QPT [U HH1 fTT] - However, their role in QPT has been 
limited to the discrete variable state space. The technique presented in this paper 
extends the purview of maximum-likelihood QPT to the continuous variable state space, 
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Figure 1: Schematic of the experimental setup for performing csQPT. 

thereby allowing physically consistent quantum process estimation through homodyne 
tomography experiments [12]. A further advantage of the present technique is the need 
of a significantly narrower range of coherent states to probe the process as compared 
to [5j [6]. We test our approach on a number of processes that are relevant to quantum 
optical information processing: identity, attenuation and photon creation. In doing so, 
we elaborate a number of recommendations for practical use of the method. 

2. The method 

2.1. Iterative process estimation using Jamiolkowski isomorphism 

The process reconstruction scheme presented in this paper is based on application of 
a maximum-likelihood based QPT scheme [H [TT] to quadrature measurements in the 
Hilbert space associated with a harmonic oscillator. Consider a quantum optical process 
£ acting upon an optical mode prepared in some quantum state p m . The positivity 
of density matrices deems it necessary that £ be a completely positive (CP) map, 
in addition to being trace non-increasing [13]. The output state £(p m ) of such a 
process can be subjected to optical homodyne measurements of its field quadratures 
xq = x cos 9+psin 9, where x and p are the canonical position and momentum operators 
and 9 is the local oscillator phase. For the output of the probe p m , the probability of 
detecting a specific quadrature value x for a phase 9 is given by 



where 11(9, x) = \9,x) (9,x\ is the projector associated with the quadrature eigenstate 
\9,x) and the superscript m on the left hand side denotes the probe state index. The 
above expression can be considered as a probability distribution function with £ as 
the parameter. If one performs N measurements for each of the M input probe states 
p m , obtained as a set of phase and corresponding quadrature values {9 itm , x^ m } where 
1 < i < N and 1 < m < M, one can obtain the log-likelihood functional as 




(1) 



= £ In (^>, m )) 
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= J] In (Tr [ft(0 iim , x i>m )f (p m )] ) . (2) 

This functional is convex over the space of CP maps [Hj. The objective of maximum- 
likelihood estimation is to determine the parameter £ est that is as close to the actual 
parameter as possible, by maximizing the likelihood functional £(£) over the space of 
CP maps 

£ est = argmax£(£ ). (3) 

This optimization problem is not straightforward and has been handled previously 
through various methods such as the uphill simplex [15]. However, a more rigorous 
yet technically simpler approach involves the formulation of an extremal equation that 
maximizes the log-likelihood functional given in equation (J2J). 

In order to carry out the reconstruction procedure, one needs to first select a certain 
basis for the representation of the process and the relevant operators. In the Fock 
(number state) basis, the quantum process can be represented by a rank-4 tensor that 
relates the density matrix of the input and output states as [5, 6J 

[Pout] jk = £jk U [Pin} mn , (4) 

m,n,j,k 

where £^ n = (j\£ (\m) (n|) \k) and p mn = (m\p\n). Although the optical Hilbert space 
is of infinite dimension, in practical process tomography it is truncated to the spanning 
set of several lowest Fock states, as will be discussed later. Also, the projectors 11(8, x) 
can be expressed in this basis as 

n mn (0,x) = (m\fl(8, x)\n) = (m\9,x) (8,x\n) , (5) 

where the overlap of the quadrature eigenstate with the number state is given by [TUl HE] 

<*-**> fey ^k* ■■ w 

With the selected basis, we proceed to formulating a numerical procedure for the 
reconstruction of the quantum process. For a concise mathematical visualization, we 
resort to the Jamiolkowski isomorphism between linear CP maps £ from operators on 
the Hilbert space H to the space /C and positive semidefmite operators E on the Hilbert 
space Ti <S> JC. The explicit relation between E and £ is given as [14] 

E= £™ l m > H ® \3) W ■ ( ? ) 

m,n,j,k 

With the definition in equation (J7]), the output p out G /C of a process £ for an input 
Pin e H is 

Pout = £{p- m ) = Tr w [EpT ® 4] , (8) 

where T denotes transposition in the number basis. In addition, one must apply the 
trace-preservation condition (Tr[p out ] = Tr[p in ]) over the process £, which yields 

Ti K [E] = I H . (9) 
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The reconstruction procedure can be extended to also encompass trace non-preserving 
processes, as will be shown subsequently. The problem has thus reduced to the 
determination of (diniHdim/C) 2 parameters subject to diniH 2 constraints. When the 
input and output Hilbert spaces are identical, this amounts to evaluating dim"H 4 — dim?/ 2 
free parameters. 

For the process output of the input probe state p m , the probability of reading a 
quadrature value x for a given local oscillator phase 6 can be obtained by substituting 
(IE]) into equation (pQ) to obtain 

E Pm 



pT(x) 



Tr 



IL(9,x) 



(10) 



Operator E should then maximize a constrained log-likelihood functional in order to 
stand as the most likely quantum process that has the set of outcomes {0 itrn ,x itm } for 
the input probes {p m }- The relevant log-likelihood functional is given as 



£(E) = ^ln(< m (x, 

m,i 

= ln ( Tr Ep 



Ti[AE] 



i.mj X-i^ m 



Tr[A£], 



11^ 



m,i 



where A = A(§>/k: and A is the Hermitian matrix of Lagrange multipliers that incorporates 
the trace preservation condition ([9]). Again, and x i)m belong to the set of quadrature 
data for the m th probe state given by {9i ym , Xj /m }. An extremal equation can be obtained 
by varying equation (TlTI) with respect to E: 

5£{E) = C(E + 5E) - £(E) = 0, (12) 

which gives 

1 



Tr 



E 



<g> U.(9 itTn , x i)m ) - A ) 5E 



0. 



(13) 



This holds for all 5E, so that the expression in the parentheses can be equated to zero 
and one has 



where 



E = A' 1 RE, 



m,i 



(•^i,m) 



Pm 



i,m.) •Ei,m 



(14) 



(15) 



Owing to Hermicity, one may rewrite equation ( TT4"|) as E = ERA 1 . Using this, along 
with equation (TT4T) . we arrive at 

E = A^RERA- 1 . (16) 

A can be determined by substituting the expression for E in equation ( Fl6i) into the 
trace-preservation condition Qj: 



A= [Tr K [RER] 



1/2 



(17) 
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Equations (1161) and (1171) can be solved numerically through iterations, starting from an 
unbiased initial E, such as E^ = /^^/(dim/C). At each step of the iterations, the 
positive semi-defmiteness of E is ensured and the constraint Tr^f-E] = In is satisfied. 

Quantum processes may also be probabilistic, in which case the trace of the input 
quantum state is not preserved. The probability of occurrence of a probabilistic quantum 
process is given by 

Psucccss = Tr [£(p)], (18) 

The reconstruction of probabilistic quantum processes can be viewed as a reconstruction 
of a trace-preserving, deterministic CP map £ if the failure of the process is taken to be 
a measurement event associated with the projection operator fl© onto a fictitious state 
|0) [14J. In order to analyze such a process, one can extend the Hilbert space to form 
^Qotai = /C © /Cf a ii, where /Cf ai i is spanned by the single state |0). The original set of 
projectors fl$(x) for each 6 is augmented by adding IT0 so that the new set of projectors 
satisfies the closure relation over /C to tai, i-e. W Jflg(x)dx + 11$ = I. Subsequently, the 
likelihood functional, with the extended trace-preserving map £ as parameter, can be 
rewritten as 

C(E) = [ftnln (iCfom)) + ( X - 9m)^ (pJT)] " Tr [ A 4 (19) 

m,i 

where g m is the fraction of successful events over total events, which can be determined 
experimentally. The extremal equation would then contain a modified operator R given 
by 



m,i 



9m ~T ^ tt f a \ i 9m ~T 



p m ® u(e i>m , x i>m ) + — —p m ® n 

■^i,m) F0 



(20) 



Iterations can now be performed with the new R to obtain the trace-preserving process 
tensor E. The actual process tensor E is obtained by taking the projection of the 
estimated tensor E onto the subspace H® K. 

Our analysis so far did not specify which states were to be used as probes; the only 
requirement is that these states compose a spanning set in the space of density matrices. 
In csQPT, the role of probe states is played by coherent states [5]. The density operator 
of an arbitrary state can be written as a linear combination of coherent state density 
operators using the optical equivalence theorem: 

p = J Pp{a) \a) (a\d 2 a, (21) 

where Pp{a) is the Glauber-Sudarshan P function of state p. Using the linearity of 
quantum processes with respect to density matrices, the process output is then given 
by 

£{p) = J P p £{\a) (a\)d 2 a. (22) 

Therefore, if the response of the quantum system to all coherent states is known, the 
output of any arbitrary unknown quantum state can be computed. In other words, 
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measurements on the set of responses £(\a m ) (a m \) for coherent states \a m ) provides 
tomographically complete information about the quantum process. 

2.2. Practical issues 

We now proceed to discussing a few practical issues arising in the implementation of 
the above algorithm of csQPT. The first issue is associated with infinite dimension of 
the optical Hilbert space. In practical implementation of csQPT, the process tensor is 
reconstructed for a subspace H(n max ) of the Hilbert space spanned by Fock states up 
to a certain cut-off value, n max . The choice of n max is correlated with the maximum 
amplitude a max of the set of coherent probe states. Given a data set with a specific 
a max , the choice of n max depends on many factors, in particular, the process itself (see 
supplementary online material to [5]). 

For the iterative cycle, the cut-off value must be chosen sufficiently high so that 
H(n max ) accommodates all of the coherent probe states and the associated output states. 
Otherwise, the probe states and the quadrature data will be inadequately represented 
by H(n max ). This will lead to inaccurate reconstruction of the process tensor; we refer 
to this phenomenon as truncation errors. 

For physically realistic processes, we expect the fractions of \a m ) and £(\a m ) (a m \) 
that lie outside the reconstruction subspace to vanish as n max tends to infinity. Hence, for 
a given a max , it is possible to choose a value of n max such that the associated truncation 
errors are arbitrarily low [17] . 

However, a high cut-off value may give rise to another class of inaccuracies, which 
we call data insufficiency errors. If the overlap of a given Fock state \n) with all of 
the \ct m ) is low, so is the contribution of |n) to the log-likelihood functional, and hence 
the available data will not provide sufficient information about the effect of the process 
on \n). In contrast to the truncation errors, the data insufficiency errors grow with 
n max , but only apply to the process tensor elements associated with high input photon 
numbers. 

Therefore the following dual-step procedure for the choice of the cut-off may be 
necessary. The initial value of n max must be sufficiently high to ensure absence of 
truncation errors. Subsequently, after the iterative cycle has been completed, we choose 
a secondary cut-off value, n' m3X < fi max , and remove all the process tensor elements 
containing indices above n^ ax . The choice of n^ ax can be determined by calculating 
the statistical errors associated with each process tensor element - similar to the error 
estimations for state tomography [131 [XH1 [XH], ED] - However, further research is required 
to determine statistical errors for QPT and establish a concrete bound for n' max . In the 
next section, we illustrate the effect of the chosen subspace dimension on the process 
reconstruction through various simulations. 

As in the case of state tomography [10], our algorithm permits automatic correction 
for optical losses and inefficient detectors in the process tensor reconstruction. In order 
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to account for non-unitary efficiency r], the projection operators are replaced by 
fl v (6, x) = Y B n+k,n(jl) B m+k,m(jtl) <m|ft(0, z) \n) x |m + &) (n + k\ , (23) 

m,n,j,k 

1 /2 

where B n+k ^ n = [("^ )^ n (l — ^)' c ] • Substituting this into equation (fT5]) and 
performing the iterations generates the original process tensor pertaining to the case 
of ideal detection. 

Many physically relevant processes are phase-invariant: applying an optical phase 
shift to the input state results in the same shift to the output. Mathematically such 
processes satisfy the following relation [6j [10] 

s[u{4>)pu\4>)\ = um{p)u\<t>) (24) 

In this case, further simplifications can be made. If the action of the process on a 
coherent state \a) is known, so is the outcome for \ae l< ^). Therefore, one needs to only 
perform measurements for input coherent states with amplitudes on the positive real 
axis. When condition (I24p is applied in the Fock basis, the elements of the process tensor 
gmn £ or a ph ase _i nvar i an t process vanish except when m — n = j — k. This condition is 
incorporated into the probability distribution f fTOj) as 

p?(x) = Tr[M(E)pl®IL(e,xj\, (25) 

where M. denotes a masking operation over E. If U m = \m) (m\ denotes a projection 
operator in the number basis, then Ai(E) can be expressed as 

m{e)= <5 m _ nj -_ fc (n m (g)n i )s(n n (8)n fc ). (26) 

m,n,j,k 

Since the trace operation is invariant under cyclic rearrangements of the operators 
and the Kronecker delta is invariant under transposition of indices, the probability 
distribution ff25|) . and consequently, the expression for the operator R in equation ffT5]) 
changes to 

tf(x)=Tt[EMfa<aTL(0,xj)] (27) 
R = Y1 n m \z , M (pl ® m, m , x hm) ) . (28) 

When the above relations are used, the elements of S"pi n , for which m — n ^ j — k, 
vanish, resulting in the incorporation of the phase invariance condition. 

In some cases, the value of the log-likelihood oscillates before converging to the 
maximum owing to overshoots. Stabilization can be achieved using the diluted algorithm 
that slows down but guarantees convergence [2T]. The operator R, in that case, is 
modified to a weighted sum of itself and the identity operator as 

R! = pR + (1 - p)I, (29) 

where < p < 1. As the value of p decreases, the algorithm becomes more and more 
dilute, resulting in increased stability but a reduced rate of convergence. In addition, 
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monotonic increase of the likelihood is guaranteed for small values of \i (see Appendix A). 



One may try to find the optimal value of \i that maximizes the increase in likelihood 
at the cost of increased computational complexity Gradually varying the value of /i 
during the iterations may be justified for some processes. 

The number of quadrature measurements for each probe state typically ranges in 
tens of thousands. With multiple probe states, the iteration cycle may require significant 
computation time. In order to speed up the computation, binning of the data points in 
the quadrature and phase axes may be useful. A suitable step size is chosen for each 
trade-off between the desired computational time and the quantization error. 
For each 8(\a m ) (a m \), quadrature data points are then clubbed into bins with centers 
(Vi V)' With this modification, the log-likelihood functional ([2]) now reads as 

= h m;n,v^ PT v , m ( x u,m) , (30) 

m,u,v 

where h m . UtV denotes the number of data points in the bin with center (9 u , m , x v>m ). 
Ideally, one must obtain the POVM associated with the bin center as a function of 
all the POVMs lying in the bin. However, given a small size of the bin, this element 
can be approximated by projection onto the quadrature value at the center of the bin. 
Similarly, the operator R in equation ( TT5|) can be rewritten as 



R= Y1 m m p V \ P™ ® n (^,m,3^,m)- (31) 



For further speedup, one can compute R in a parallel fashion on different threads owing 
to absence of interdependency in the summation procedure. 

In practical experiments on probabilistic processes, the frequency of successful 
events can be low: g m <C 1. In this case, the process tensor elements of interest (i.e. 
those related to /C) will be small and thus suffer from increased relative error. This issue 
can be resolved by rescaling the values of all g m by the same factor for all probe states, 
keeping in mind the requirement that g m < 1 for all m. Physically relevant elements of 
the process tensor will then rescale by the same factor, reducing the relative error. 



3. Implementation and results 

In order to test the algorithm, we have implemented it using Matlab and studied the 
reconstruction of a few quantum processes using simulated data. Theoretical process 
tensors of identity, attenuation and photon creation [6] were used to find the marginal 
probability distribution functions for various probe states using equation (fit)]) . From the 
marginals, we generated synthetic experimental data through Monte-Carlo simulations. 

Each process was applied to four coherent probe states with a m ranging from to 
0.9375 in steps of 0.3125. For each input probe state, the output state dataset consisted 
of 100,000 phase and quadrature points {9,x}. This set of data was subjected to the 
reconstruction method described. The iterations were halted when the change in process 
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(d) (e) (f) 

Figure 2: Comparison of the diagonal values (^ m ) of theoretical and reconstructed 
process tensors, (a), (b) and (c) show process tensors for theoretical identity, attenuation 
(by factor 0.9) and photon creation processes, respectively, (d), (e) and (f) show 
the process tensors for the corresponding reconstructed processes using the algorithm 
presented in this article. The photon creation process tensor has been scaled to match 
the theoretical one (i.e. g — 1). 

tensor elements was insignificant over a large number of iterations. However, a better 
approach would be to set a threshold for the increase of the log-likelihood [22J. 

The result obtained by running the reconstruction technique is a 4-dimensional 
process tensor whose diagonal elements have a simple interpretation. For a given 
quantum process S, the diagonal element S™™ denotes the probability that the output 
contains k photons when the process is subjected to m input photons. A comparison 
between the diagonal elements of the theoretical and reconstructed process tensor is 
made in figure [2] and exhibits close match between the two. 

The process of photon creation aJ requires additional discussion because it 
corresponds to a non-unitary, trace non-preserving operator. Therefore, in experimental 
practice it can only be implemented probabilistically. The optical mode containing the 
target state \tp) is directed into the signal channel of a parametric down-conversion setup. 
The state of the down-conversion output in the signal (s) and idler (i) channels can 
then be written as \ip) s |0) i + g(a) \ip)) s where g is the down-conversion amplitude. 
Detection of a photon in the idler channel projects the signal state onto a) thereby 
heralding a photon addition event [23]. For coherent state input \ip) = \a), the 
event probability, corresponding to the quantity g m in equation (fl9l) . is proportional 
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(a) (b) 

Figure 3: Reconstruction of the photon creation process with correction for inefficiency: 
(a) r\ = 0.75 and (b) rj = 0.55. 

to \g\ 2 (l + \a\ 2 ). 

We take the value g 2 = 0.1 during simulations to ensure that success probabilities 
remain less than 1 for the probe states selected. This makes the photon creation 
process trace non-increasing and thus physical. Note that the process tensor reported 
in figure Mf) has been normalized by dividing by g 2 , so that its scale matches that of 
the process tensor for the photon creation operator a) given in figure [2]^c) . 

The iterative reconstruction of photon creation exhibited relatively poor 
convergence. Diluted iterations (|29|) were required in the beginning in order to curb 
oscillations. However, as the iterations progressed, the rate of increase of the likelihood 
value became extremely low. We circumvented this issue by implementing the successive 
over-relaxation technique. Setting fi in equation ( 129]) to slightly over 1 while iterating 
accelerated the increase in likelihood. As soon as a decrease in the likelihood value 
was registered due to an overshoot, ji was reset to 1, and then slowly increased again 
after stabilization. This procedure was applied multiple times until a fair amount of 
convergence was observed. In a loose sense, the over-relaxation method employs linear 
extrapolation by selecting a tensor that lies on the line joining the current iterate and the 
next iterate but is beyond the latter by a fraction. If the iterations happen to proceed in 
the direction of maximum likelihood gradient, it allows faster convergence by inducing 
greater leaps. Additionally, it may also help in escaping limit cycles encountered during 
the iterations. 

We have also tested the reconstruction technique for photon creation in the case of 
inefficient detection. The output density matrices for the probe states were calculated 
using the beam splitter model of absorption [16] . With these modified density matrices, 
we have generated test data using Monte Carlo simulations for rj = 0.75 and rj = 0.55. 
A comparison of the reconstructed process tensors is given in figure [3] 

Finally, we investigated the effect of the dimension of the subspace of optical Hilbert 
space chosen for the reconstruction, specifically for the photon creation process. In 
order to eliminate statistical errors in this reconstruction, we directly used the marginal 
distributions instead of simulated quadrature data sets to obtain the values of h m - UiV in 
equations ( 130]) and (|3T]) . The performance criterion is taken to be the worst-case fidelity 
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(a) 

Figure 4: Effect of the photon number cut-off (n max ) on the reconstruction of the 
photon creation process, (a) Worst-case process reconstruction fidelity as a function 
of n max for a max = 0.4 (top) and a max = 0.6 (bottom) and different n' m3X . The slight 
decreases of fidelity with increasing n max and constant n' max , observed in some cases, are 
numerical artefacts, (b) Diagonal values of the process tensor S™™ for a max = 0.6 and 
^max = 3. The reconstructed process tensor has significant artefacts due to truncation 
errors, (c) Diagonal values of the process tensor £j^ n for a max = 0.6 and n max = 8. The 
reconstructed process tensor elements associated with input photon numbers 6, 7 and 8 
are invalid due to data insufficiency (for example, |(a m ax = 0.6|n = 6)| 2 = 2.1 x 10~ 6 ). 



[24] over the input space H(n max — 1), defined as 

F(£,£ cst ) = . min Tr U y/£{fi £ cst (p) y/e(fi) , (32) 

pGH(n max -l) y ' J 

where £ is the actual process tensor and £ est is the estimated process tensor. Note that 
the photon number cutoff for the fidelity calculation is taken to be n max — 1 to ensure the 
Hilbert space closure under photon addition. Minimization over EI(n max — 1) is carried 
out through a Monte Carlo simulation that involves introducing small random changes 
in the density matrix p within EI(n max — 1) and accepting the change whenever the value 
of the fidelity decreases. 

The solid line in figure 11(a) shows the worst-case fidelity versus n max for two values 
of a max . For each given a max , the fidelity initially increases with n max as the truncation 
effects subside and decreases afterwards due to data insufficiency The range of n max , 
over which the process tensor is reconstructed correctly, shifts towards the higher photon 
numbers with increasing a max owing to greater contribution of higher photon numbers 
in probe states of higher amplitudes. 

Figure BJb,c) further illustrates the two types of errors associated with the choice 
of the cut-off point. If n max is chosen too low (figure BJb)), truncation errors 
compromise the entire reconstructed process. If the reconstruction subspace is sufficient 
to accommodate all the input probe states and associated output states (figure 11(c)), 
only the process tensor elements associated with high input photon numbers are 



Maximum-likelihood coherent- state quantum process tomography 



13 



reconstructed incorrectly. In this case, introducing a secondary cut-off at n' mSuX = 5 
is justified. 

The dashed lines in figure HJ^a) display the advantages of the dual cut-off approach 
introduced in section 12.21 For example, with a max = 0.6, optimal reconstruction is 
attained with the initial cut-off point n max > 6 and subsequent cropping of the process 
tensor with n' max = 5. With this approach, the worst-case reconstruction fidelity is 
higher than for all cases with n' max = n max . Note, however, that in most examples we 
studied, the dual cut-off method offers only a small advantage and may not be justified 
in practical csQPT. 

As evidenced by figure Eta), the secondary cut-off points should be chosen close to 
n max = 3 and 5, for a max = 0.4 and 0.6, respectively. These values are much higher than 
those calculated in the supplementary information to Ref. [5]. Specifically, the method 
of Ref. [5] for the same values of n max = 3,5 would require the maximum coherent state 
amplitudes of 8 and 12, respectively. Further, the method of Ref. [6] requires 2N probe 
states for reconstruction in a Fock space of dimension d = N+l, while our method poses 
no such constraints. In other words, for a given set of probe coherent states (defined 
by their number and maximum amplitude a max ), the present reconstruction method 
provides much more information about the process tensor than previous methods. 

One must note that further research is needed for the inverse problem of determining 
the optimum ct max for a chosen Fock space dimension d = n max . According to the 
numerical examples we studied, it is reasonable to choose a max such that (a max |n max ) 2 « 
1/N where N is the total number of quadrature measurements. 

4. Summary 

We have presented a maximum-likelihood based experimental data processing technique 
for the tomographic reconstruction of quantum optical processes. This technique relies 
on measuring the response of the process to various coherent probe states through 
optical homodyne tomography. The reconstruction applies directly to the obtained 
data, unlike the previous coherent state QPT methods that involve intermediate 
reconstruction of density matrices of the output states. The range of probe states 
required for reconstruction has also been reduced. Complete positiveness and trace 
preservation/non- increase conditions are incorporated in the estimated process tensor 
by imposing a priori constraints, thus yielding physical results. The simplicity and 
robustness of this technique make it appealing for quantum process estimation, with 
applications extending to optical quantum computing and quantum communication. 
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Appendix A. Proof of monotonic increase of Log-likelihood for the diluted 
algorithm 



In this section, we shall prove that the diluted algorithm of equation f[29|) ensures 
monotonic increase of the log-likelihood value for < \x << 1. We start by considering 
the (k + l) th iteration 



E {k+1) - A { ^R {k) E {k) R {k) A^y 
As per the diluted algorithm, R^ is modified as 

R'(k) — ^R{k) + (1 _ A i )^cg)/C- 
To find the an expression for the normalization operator A', 



X', 



(k) 



(A.l) 

(A.2) 
In = 



(Tryc[-R( fc) -E(fc)-R( fc) ]) 1/2 <8> In, we first evaluate R' (k) E( k )R' (k) to first order in /i: 



R{k) E (k)R{k) 



(ky 



\k) 



fiRw + (i-fi)i 



n®K 



E 



n®K 



E(k) — 2fiE^ + fj,R( k )E(k) + fiE(k)R\k) + 0(fi 2 ). (A. 3) 



The matrix X', u \ can be obtained as 



a; 



(fc) 



(Tr K [R[ k) E {k) R[ k) ]y/ 2 



Tr/c 



E( k ) — 2[iE(k) + /i lR(k)E(k) + i?( fc )-R(fc) 



1/2 



;i - 2fj.)I H + 2/iTr^ 



R(k)E( k ) + E( k )R 



-1 1/2 



(k)^(k) 



(A.4) 



(A 



' \-i 



1 - + /j,Tr K 



1 + //)/ w - /xTr^ 



R(k)E( k ) + E(tyR(ty 




(A.5) 



where in equation (1A.4I) . we have used the trace preserving condition from equation (j9]) 
Thus, one has 



(A'r 



where Yr 



(*) 



(l + li)In®K - V>Y(k), 

R(k)E( k ) + 



(A.6) 



Tr 



/c 



Equation flA.lj) can now be written as 



( A + I i "= ( A '(fc) ) 1 ^(fe) R [k) (^(fc) ) 



[(1 + h)Iu®k — vY(k)][E{k) — 2fxE^ + fj,R( k )E( k ) + fj,E( k )R^} 
E (k) + AE (k) , (A.7) 



Maximum-likelihood coherent- state quantum process tomography 15 

where AE (k) = fj, ^R {k) E {k) + E {k) R (k) - Y {k) E {k) - E {k) Y {k ^j. The log-likelihood at the 
(k + l) th iteration is given by 

C(E (k+1) ) =J2 ln ( Tl [(A*) + Mfc))^™ ® n(Zi, m AJ 

Tr 



+ J> 1 + 



m,i 



Tr 



= J^ln^Tr %)^®n(i !im ,^ m ) ) 



m,i 



Tr 



Tr 



£(fc)/5£ ® II(a; i ,mi v i,m) 



C(E {k) ) + Tr[AE (k) R {k) } 



(At 



To prove the monotonicity of the log-likelihood functional, we require that 



Tr 



AE( k \R 



(fc)-n-(fc) 



iteration, i.e. 



^Tr R( k )E( k )(R( k ) - F( fe )) + (_R (fc) - Y^)E^R^ 



> holds at each 



Tr 



R(k)E(k)(R(k) - Y(k)) + (R(k) - Y( k ))E( k )R( k ) > V k. 



(A.9) 



Starting with the left hand side (dropping subscripts): 



Tr 



RE(R — Y) + (R — Y)ER 



= Tr 
= Tr 



2RER - REY - YER 



2RER - 2REY - 2YER + (REY + YER)] . (A. 10) 



Considering the expression (REY + YER): 



Tr 



REY + YER 



= Tr 
= Tr 

= Tr 
= Tr 



(RE + ER)Y 
n RE + ER 



Tr 



/ RE + ER 



IK 



2 Tr, 



2 Tr 



f RE + ER 




Tr K (E) 



(A.ll) 
(A.12) 
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K. 



(A.13) 



where in flA.llj) . we have used a property of the partial trace Tr^: Tr A(B ® 1, 



</cJ 



Tr 



BTi/c(A) and in flA.12j) . we have used equation ([9]). Substituting (IA.13j) in equation 



f lA.lOp . we have 



Tr 



RE(R -Y) + (R- Y)RE 
= 2Tr \RER - REY - YER Y EY 
= 2Tt[(rE 1 / 2 -YE 2 ) {E V2 B- E' 2 Y 

= 2Tr 



(A. 14) 
(A.15) 

= 2Tr [X^X] > where X = E 1/2 R - E 1/2 Y, (A. 16) 

where in (1A.14|) . the positive semidefiniteness of E allows us to factorize it as E — 



(e 1 ' 2 R - E^fV (e 1 ' 2 R - E^Y 



E X I 2 E X I 2 ! , with E l l 2 being a Hermitian matrix. ( 1A.15j) follows as R and Y are also 
Hermitian matrices. We arrive at the inequality (IA.16I) as the matrix X^X is positive 
semidefinite and thus has non-negative trace. This completes the proof that the log- 
likelihood monotonously increases with k for < [l « 1. 
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